Author

Victor L. Jardim

System and models

Real system

Two-Way Econometric Fixed Effects model

Group Mean Covariate Model # Richness models

Naive, mixed effects and fixed effects models

Code
tar_load(alphamod)
alphamod <- alphamod %>% mutate(across(c(PC1_score, PC2_score, Depth, Fetch_max, Current_mean, T_mean, OM), scale)) %>% mutate(SY = paste(.$Site, .$Year, sep = "_"))


mod_naive <- lm(S.obs ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean, data = alphamod)

mod_re <- lmer(S.obs ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+(1|Site), data = alphamod)

mod_fe <- lm(S.obs ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+Site+Year, data = alphamod)

mod_twfe <- lm(S.obs ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+SY, data = alphamod)
Code
export_summs(
  mod_naive,
  mod_re,
  mod_fe,
  mod_twfe,
  model.names = c("naive", "re", "fe", "twfe" ),
  scale = TRUE,
  statistics = c(
    N = "nobs",
    DF = "df.residual",
    AIC = "AIC",
    BIC = "BIC",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
  coefs = c(
    "PC1" = "PC1_score",
    "PC2" = "PC2_score",
    # "PC1_site" = "PC1_site",
    # "PC2_site" = "PC2_site",
    # "PC1_sy" = "PC1_sy",
    # "PC2_sy" = "PC2_sy",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Fetch" = "Fetch_max",
    "Current" = "Current_mean",
    "Organic matter" = "OM",
    "Temperature" = "T_mean"),
  robust =  "HC1", error_pos = "right"
)
naiverefetwfe
PC12.46 *  (1.13)1.63   (1.33)-0.12   (1.15)0.20    (0.88)
PC24.26 ***(1.14)3.82 * (1.63)2.71   (1.45)2.08    (1.19)
PC1:PC20.43    (0.90)2.30 * (1.11)2.67 * (1.04)2.83 ***(0.82)
Fetch-2.52    (1.39)0.25   (3.27)13.90 **(4.39)13.33 ***(3.86)
Current-1.10    (1.16)-4.60   (3.45)-25.46 **(8.62)-18.04    (10.89)
Organic matter2.97 ** (0.95)4.07 **(1.23)2.30 * (0.98)3.62 ***(1.08)
Temperature2.44 ** (0.83)1.81 * (0.84)-1.81   (2.07)-16.84    (45.33)
N348           348          348          348           
DF339.00        337.00       319.00       222.00        
AIC2908.35        2873.04       2785.52       2712.69        
BIC2946.88        2915.42       2901.09       3201.92        
Marginal R2           0.26                            
Conditional R20.43        0.45       0.64       0.83        
Adjusted R20.42                  0.61       0.74        
F32.04                  20.55       8.94        
p-value0.00                  0.00       0.00        
All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. Standard errors are heteroskedasticity robust. *** p < 0.001; ** p < 0.01; * p < 0.05.
Code
(naivecoef <- coeftest(mod_naive, 
         vcov = vcov(mod_fe,
                       # cluster =  ~ Site,
                       type = "HC1")) |>
tidy() |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "naive"))
# A tibble: 7 × 6
  term                estimate std.error statistic p.value model
  <chr>                  <dbl>     <dbl>     <dbl>   <dbl> <chr>
1 PC1_score              2.46       1.21     2.03  0.0432  naive
2 PC2_score              4.26       1.51     2.82  0.00513 naive
3 OM                     2.97       1.21     2.46  0.0145  naive
4 T_mean                 2.44       2.17     1.13  0.261   naive
5 Fetch_max             -2.52       4.80    -0.525 0.600   naive
6 Current_mean          -1.10       9.65    -0.114 0.910   naive
7 PC1_score:PC2_score    0.425      1.04     0.408 0.683   naive
Code
(recoef <- tidy(mod_re) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "naive")) %>% 
  select(-effect:-group)
# A tibble: 7 × 5
  term                estimate std.error statistic model
  <chr>                  <dbl>     <dbl>     <dbl> <chr>
1 PC1_score              1.63      1.32     1.24   naive
2 PC2_score              3.82      1.57     2.44   naive
3 OM                     4.07      1.21     3.35   naive
4 T_mean                 1.81      0.837    2.17   naive
5 Fetch_max              0.249     3.08     0.0810 naive
6 Current_mean          -4.60      3.34    -1.38   naive
7 PC1_score:PC2_score    2.30      1.10     2.09   naive
Code
coeftest(mod_fe, 
         vcov = vcovCL(mod_fe,
                       cluster =  ~ Site,
                       type = "HC1")) |>
  tidy() |>
  filter(term %in% c(
    "PC1_score",
    "PC2_score",
    "PC1_score:PC2_score",
    "Fetch_max",
    "Current_mean",
    "OM",
    "T_mean"))
# A tibble: 7 × 5
  term                estimate std.error statistic  p.value
  <chr>                  <dbl>     <dbl>     <dbl>    <dbl>
1 PC1_score             -0.118     1.10     -0.107 9.15e- 1
2 PC2_score              2.71      1.12      2.42  1.62e- 2
3 OM                     2.30      0.966     2.39  1.76e- 2
4 T_mean                -1.81      3.65     -0.496 6.20e- 1
5 Fetch_max             13.9       2.07      6.73  7.98e-11
6 Current_mean         -25.5       9.17     -2.78  5.83e- 3
7 PC1_score:PC2_score    2.67      1.09      2.44  1.51e- 2
Code
(fecoef <- coeftest(mod_fe, 
         vcov = vcovCL(mod_fe,
                       cluster =  ~ Site + Year,
                       type = "HC1")) |>
  tidy() |>
  filter(term %in% c(
    "PC1_score",
    "PC2_score",
    "PC1_score:PC2_score",
    "Fetch_max",
    "Current_mean",
    "OM",
    "T_mean")) |>
  mutate(model = "fe"))
# A tibble: 7 × 6
  term                estimate std.error statistic    p.value model
  <chr>                  <dbl>     <dbl>     <dbl>      <dbl> <chr>
1 PC1_score             -0.118     1.12     -0.106 0.916      fe   
2 PC2_score              2.71      1.02      2.66  0.00810    fe   
3 OM                     2.30      0.975     2.36  0.0187     fe   
4 T_mean                -1.81      3.06     -0.592 0.555      fe   
5 Fetch_max             13.9       2.94      4.73  0.00000334 fe   
6 Current_mean         -25.5      10.1      -2.52  0.0122     fe   
7 PC1_score:PC2_score    2.67      1.02      2.61  0.00944    fe   
Code
(twfecoef <- coeftest(mod_twfe, 
         vcov = vcovCL(mod_twfe,
                       cluster =  ~ Site + Year,
                       type = "HC1")) |>
tidy() |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "twfe"))
# A tibble: 7 × 6
  term                estimate std.error statistic   p.value model
  <chr>                  <dbl>     <dbl>     <dbl>     <dbl> <chr>
1 PC1_score              0.202     1.29      0.156 0.876     twfe 
2 PC2_score              2.08      1.49      1.40  0.163     twfe 
3 OM                     3.62      0.824     4.39  0.0000177 twfe 
4 T_mean               -16.8      67.8      -0.248 0.804     twfe 
5 Fetch_max             13.3       3.23      4.13  0.0000522 twfe 
6 Current_mean         -18.0       8.96     -2.01  0.0453    twfe 
7 PC1_score:PC2_score    2.83      1.27      2.23  0.0269    twfe 
Code
coefs <- naivecoef %>% bind_rows(recoef, fecoef, twfecoef) %>% 
  mutate(term = factor(term, levels = c("PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")),
      model = factor(model, levels = c("naive", "re", "fe", "twfe")))

ggplot(coefs, aes(x = estimate, y = term, colour = model, shape = model))+ 
  geom_pointrange(aes(xmin = estimate-std.error, xmax = estimate+std.error), position = position_dodge2(.5, reverse = TRUE))+
  scale_y_discrete(limits = rev)+
  theme_light()

Naive, Grouped Mean Covariate, and Two-way grouped mean covariate models

Code
alphamod <- alphamod %>%
  group_by(Site) %>% 
  mutate(PC1_site = mean(PC1_score),
         PC2_site = mean(PC2_score)
         ) %>% 
  ungroup() %>% 
  group_by(Year) %>% 
  mutate(PC1_y = mean(PC1_score),
         PC2_y = mean(PC2_score)) %>% 
  ungroup() %>% 
  group_by(Point) %>% 
  mutate(PC1_point = mean(PC1_score),
         PC2_point = mean(PC2_score)) %>% 
  ungroup()
  
mod_gmcov <- lmer(S.obs ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + Depth+OM+T_mean+Fetch_max + Current_mean+ (1|Site) + (1|Year), data = alphamod)
mod_fegmcov <- lmer(S.obs ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + Depth+OM+T_mean+Fetch_max + Current_mean+ (1|Site) + Year, data = alphamod)
mod_twgmcov <- lmer(S.obs ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + PC1_y + PC2_y + Depth+OM+T_mean+Fetch_max + Current_mean + (1|Site) + (1|Year), data = alphamod)

export_summs(
  mod_naive,
  mod_gmcov,
  mod_fegmcov,
  mod_twgmcov,
  model.names = c("naive", "gmcov", "fegmcov", "twgmcov" ),
  scale = TRUE,
  statistics = c(
    N = "nobs",
    DF = "df.residual",
    AIC = "AIC",
    BIC = "BIC",
    `Marginal R2` = "r.squared.fixed",
    `Conditional R2` = "r.squared",
    `Adjusted R2` = "adj.r.squared",
    `F` = "statistic",
    `p-value` = "p.value"
  ),
  coefs = c(
    "PC1" = "PC1_score",
    "PC2" = "PC2_score",
    "PC1_site" = "PC1_site",
    "PC2_site" = "PC2_site",
    "PC1_y" = "PC1_y",
    "PC2_y" = "PC2_y",
    "PC1:PC2" = "PC1_score:PC2_score",
    "Fetch" = "Fetch_max",
    "Current" = "Current_mean",
    "Organic matter" = "OM",
    "Temperature" = "T_mean")
)
naivegmcovfegmcovtwgmcov
PC12.46    0.43   0.41   0.42   
(1.27)   (1.21)  (1.21)  (1.21)  
PC24.26 ***2.43   2.44   2.44   
(1.16)   (1.52)  (1.51)  (1.52)  
PC1_site       13.14 * 12.77 * 13.16 * 
       (5.65)  (5.78)  (5.67)  
PC2_site       11.63 * 11.69 * 11.62 * 
       (4.81)  (4.93)  (4.82)  
PC1_y                   2.70   
                   (3.53)  
PC2_y                   -0.88   
                   (3.64)  
PC1:PC20.43    3.02 **3.02 **3.02 **
(0.98)   (1.00)  (1.00)  (1.00)  
Fetch-2.52    10.60 * 10.82 * 10.63 * 
(1.43)   (4.43)  (4.46)  (4.43)  
Current-1.10    -4.76   -5.65   -4.87   
(1.15)   (4.58)  (4.74)  (4.61)  
Organic matter2.97 ** 3.30 **3.16 **3.19 **
(1.11)   (1.16)  (1.17)  (1.17)  
Temperature2.44 ** 0.01   -1.75   -0.10   
(0.85)   (1.72)  (2.17)  (1.81)  
N348       348      348      348      
DF339.00    334.00   324.00   332.00   
AIC2908.35    2799.41   2739.44   2794.70   
BIC2946.88    2853.34   2831.90   2856.33   
Marginal R2       0.36   0.46   0.37   
Conditional R20.43    0.70   0.71   0.71   
Adjusted R20.42                      
F32.04                      
p-value0.00                      
All continuous predictors are mean-centered and scaled by 1 standard deviation. The outcome variable is in its original units. *** p < 0.001; ** p < 0.01; * p < 0.05.
Code
gmcoef <- tidy(mod_gmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "gmcov") %>% 
  select(-effect:-group)

fegmcoef <- tidy(mod_fegmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "fegmcov") %>% 
  select(-effect:-group)

twgmcoef <- tidy(mod_twgmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "twgmcov") %>% 
  select(-effect:-group)


coefs <- naivecoef %>% bind_rows(recoef, fecoef, twfecoef, gmcoef, fegmcoef, twgmcoef) %>% 
  mutate(term = factor(term, levels = c("PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")),
      model = factor(model, levels = c("naive", "re", "fe", "twfe", "gmcov", "fegmcov", "twgmcov")))
ggplot(coefs, aes(x = estimate, y = term, colour = model, shape = model))+ 
    geom_vline(xintercept = 0, linetype="dashed")+
  geom_pointrange(aes(xmin = estimate-std.error, xmax = estimate+std.error), position = position_dodge2(.5, reverse = TRUE))+
  scale_y_discrete(limits = rev)+
  theme_light()

Density models

Code
dens_naive <- lm(log(Fauna_Density) ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean, data = alphamod)

dens_re <- lmer(log(Fauna_Density) ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+(1|Site), data = alphamod)

dens_fe <- lm(log(Fauna_Density) ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+Site+Year, data = alphamod)

dens_twfe <- lm(log(Fauna_Density) ~ PC1_score+PC2_score+PC1_score:PC2_score+Depth+OM+T_mean+Fetch_max+Current_mean+SY, data = alphamod)

dens_gmcov <- lmer(log(Fauna_Density) ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + Depth+OM+T_mean+Fetch_max + Current_mean+ (1|Site) + (1|Year), data = alphamod)

dens_fegmcov <- lmer(log(Fauna_Density) ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + Depth+OM+T_mean+Fetch_max + Current_mean+ (1|Site) + Year, data = alphamod)

dens_twgmcov <- lmer(log(Fauna_Density) ~ PC1_score+ PC1_site + PC1_score:PC2_score+ PC2_score+ PC2_site + PC1_y + PC2_y + Depth+OM+T_mean+Fetch_max + Current_mean + (1|Site) + (1|Year), data = alphamod)

(naivecoef <- coeftest(dens_naive, 
         vcov = vcov(dens_fe,
                       # cluster =  ~ Site,
                       type = "HC1")) |>
tidy() |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "naive"))
# A tibble: 7 × 6
  term                estimate std.error statistic    p.value model
  <chr>                  <dbl>     <dbl>     <dbl>      <dbl> <chr>
1 PC1_score             0.134     0.0493     2.71  0.00706    naive
2 PC2_score            -0.0521    0.0614    -0.849 0.397      naive
3 OM                    0.196     0.0490     4.01  0.0000760  naive
4 T_mean                0.136     0.0879     1.55  0.123      naive
5 Fetch_max             0.228     0.195      1.17  0.244      naive
6 Current_mean         -0.305     0.392     -0.778 0.437      naive
7 PC1_score:PC2_score  -0.200     0.0423    -4.74  0.00000314 naive
Code
(recoef <- tidy(dens_re) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "re") %>% 
  select(-effect:-group))
# A tibble: 7 × 5
  term                estimate std.error statistic model
  <chr>                  <dbl>     <dbl>     <dbl> <chr>
1 PC1_score             0.0541    0.0527     1.03  re   
2 PC2_score             0.125     0.0652     1.92  re   
3 OM                    0.151     0.0485     3.11  re   
4 T_mean                0.0735    0.0325     2.26  re   
5 Fetch_max             0.148     0.166      0.896 re   
6 Current_mean         -0.436     0.224     -1.95  re   
7 PC1_score:PC2_score   0.0360    0.0437     0.822 re   
Code
gmcoef <- tidy(dens_gmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "gmcov") %>% 
  select(-effect:-group)

fegmcoef <- tidy(dens_fegmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "fegmcov") %>% 
  select(-effect:-group)

twgmcoef <- tidy(dens_twgmcov) |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")) |>
  mutate(model = "twgmcov") %>% 
  select(-effect:-group)

(fecoef <- coeftest(dens_fe, 
         vcov = vcovCL(dens_fe,
                       cluster =  ~ Site + Year,
                       type = "HC1")) |>
  tidy() |>
  filter(term %in% c(
    "PC1_score",
    "PC2_score",
    "PC1_score:PC2_score",
    "Fetch_max",
    "Current_mean",
    "OM",
    "T_mean")) |>
  mutate(model = "fe"))
Warning in sqrt(diag(se)): NaNs produced
# A tibble: 7 × 6
  term                estimate std.error statistic p.value model
  <chr>                  <dbl>     <dbl>     <dbl>   <dbl> <chr>
1 PC1_score             0.0295    0.0460     0.640  0.523  fe   
2 PC2_score             0.134     0.0592     2.26   0.0243 fe   
3 OM                    0.0864    0.0648     1.33   0.184  fe   
4 T_mean               -0.0869    0.124     -0.703  0.483  fe   
5 Fetch_max             0.246     0.257      0.957  0.339  fe   
6 Current_mean         -1.06      0.486     -2.18   0.0300 fe   
7 PC1_score:PC2_score   0.0309    0.0366     0.846  0.398  fe   
Code
(twfecoef <- coeftest(dens_twfe, 
         vcov = vcovCL(dens_twfe,
                       cluster =  ~ Site + Year,
                       type = "HC1")) |>
tidy() |>
      filter(term %in% c(
      "PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean"))
  |>
  mutate(model = "twfe"))
# A tibble: 7 × 6
  term                estimate std.error statistic p.value model
  <chr>                  <dbl>     <dbl>     <dbl>   <dbl> <chr>
1 PC1_score             0.0526    0.0472     1.12   0.266  twfe 
2 PC2_score             0.0954    0.0709     1.35   0.179  twfe 
3 OM                    0.139     0.0650     2.14   0.0337 twfe 
4 T_mean               -2.25      2.49      -0.905  0.367  twfe 
5 Fetch_max             0.190     0.330      0.576  0.565  twfe 
6 Current_mean         -0.596     0.592     -1.01   0.315  twfe 
7 PC1_score:PC2_score   0.0428    0.0384     1.12   0.266  twfe 
Code
coefdens <- naivecoef %>% bind_rows(recoef, fecoef, twfecoef, gmcoef, fegmcoef, twgmcoef) %>% 
  mutate(term = factor(term, levels = c("PC1_score",
      "PC2_score",
      "PC1_score:PC2_score",
      "PC1_site",
      "PC2_site",
      "Fetch_max",
      "Current_mean",
      "OM",
      "T_mean")),
      model = factor(model, levels = c("naive", "re", "fe", "twfe", "gmcov", "fegmcov", "twgmcov")))

ggplot(coefdens, aes(x = estimate, y = term, colour = model))+
  geom_vline(xintercept = 0, linetype="dashed")+
  geom_pointrange(aes(xmin = estimate-std.error, xmax = estimate+std.error), position = position_dodge2(.6, reverse = TRUE))+
  scale_y_discrete(limits = rev)+
  theme_light()

RDA

Code
tar_load(envcomp)
tar_load(bcdens)
tar_load(pal)

envcomp <- envcomp %>% mutate(Year = as.factor(Year))

list.hp4 <- list(Depth = envcomp %>% select(Depth),
                 Hydrodynamics = envcomp %>% select(Current_mean, Fetch_max),
                 Granulometry = envcomp %>% select(OM),
                 Temperature = envcomp %>% select(T_mean),
                 Complexity = envcomp %>% select(PC1_score, PC2_score),
                 Site = envcomp %>%select(Site),
                 Year = envcomp %>% select(Year))


rdamaerl.hp4 <- rdacca.hp::rdacca.hp(bcdens, list.hp4, var.part = TRUE)
upset_vp_victor(rdamaerl.hp4, cutoff = 0.005, int.var = "Complexity", pal = pal[1:8], title.cex = 22, axis.cex = 20, pch.size = 6, col.width = .8, effect.cex = 6)
Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Code
tar_read(upsetmain)

Code
tar_read(triplotrda)